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Abstract 



Numerical simulations of numerous quantum systems suffer from the 
notorious sign problem. Meron-cluster algorithms lead to an efficient 
solution of sign problems for both fermionic and bosonic models. Here 
we apply the meron concept to quantum spin systems in an arbitrary 
external magnetic field, in which case standard cluster algorithms fail. 
As an example, we simulate antiferromagnetic quantum spin ladders in a 
uniform external magnetic field that competes with the spin-spin inter- 
action. The numerical results are in agreement with analytic predictions 
for the magnetization as a function of the external field. 
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1 Introduction 



Numerical simulations of quantum systems can suffer from the notorious sign prob- 
lem. Well-known examples of systems where this problem has hindered progress are 
fermions in more than one spatial dimension, quantum antiferromagnets on non- 
bipartite lattices, as well as some quantum spin systems in an external magnetic 
field. For these systems the Boltzmann factor of a configuration in the path integral 
can be negative and hence cannot be interpreted as a probability. When the sign 
of the Boltzmann factor associated with the configuration is incorporated in mea- 
sured observables, the fluctuations in the sign give rise to dramatic cancellations. 
In particular, for large systems at low temperatures this leads to relative statistical 
errors that are exponentially large in both the volume and the inverse temperature. 
As a consequence, it is impossible in practice to study such systems with standard 
numerical methods. 

Recently, some severe fermion sign problems have been solved with a meron- 
cluster algorithm |L[ ||, [|, HI- Such an algorithm was originally developed to solve 
the sign problem associated with the phase exp(i6Q) of the Boltzmann weight in the 
(1 + l)-d 0(3) symmetric quantum field theory at vacuum angle 9 = n ||. Here Q 
is the topological charge that measures the winding number of a configuration and 
that is carried by the classical instanton solutions. The solution of the sign problem 
is based on the fact that the flip of certain clusters changes the topological charge 
by one and hence changes the sign (— 1)^ of the configuration. Such clusters can be 
associated with half a unit of topological charge and were referred to as merons Q. 
Thus meron clusters are used to identify two configurations with the same weight 
but opposite signs. Remarkably, this property of a meron is applicable to clusters 
in a variety of models with sign problems. Hence we generally use the term meron 
to denote clusters whose flip changes the sign of a configuration without changing 
its weight. In this paper we extend the meron concept to quantum spin systems in 
a magnetic field. 

As an example of physical interest we consider antiferromagnetic quantum spin 
ladders in an external uniform magnetic field that competes with the spin-spin inter- 
action. The competing interactions lead to severe problems in numerical simulations. 
The most efficient algorithms for simulating quantum spin systems in the absence of 
an external magnetic field are cluster algorithms || , in particular the loop algorithm 
|7], |8], U which can be implemented directly in the Euclidean time continuum flQ |. 



The idea behind this algorithm is to decompose a spin configuration into clusters 
which can be flipped independently with probability 1/2. This procedure leads to 
large collective moves through configuration space and eliminates critical slowing 
down. In the presence of a magnetic field pointing along the quantization axis of 
the spins, the Z(2) symmetry that allows the clusters to be flipped with probability 

1 The term meron is traditionally used to denote a half-instanton. 
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1/2 is explicitly broken. As a consequence, at low temperatures or for large values of 
the magnetic field some clusters can only be flipped with very small probability and 
the algorithm becomes inefficient. An interesting alternative to the loop algorithm 



is the so-called worm algorithm [1 ij which has produced high-precision results for 
quantum spin chains in a magnetic field [|12|]. The worm algorithm explores an en- 
larged configuration space and, unlike the loop algorithm, seems not to suffer from 
exponential slowing down at least for 1-d spin chains in a magnetic field. We are 
not aware of any applications of the worm algorithm to higher-dimensional spin sys- 
tems, and it would be interesting to see if this algorithm then still works efficiently. 
Another method based on stochastic series expansion techniques has been applied 



successfully to 2-d quantum spin systems in a magnetic field |I3| . An alternative 



strategy using cluster algorithms is to choose the spin quantization axis perpendic- 
ular to the direction of the magnetic field. Then all clusters can still be flipped with 
probability 1/2. This procedure indeed leads to an efficient algorithm for ferromag- 
nets in an external uniform magnetic field Q. Unfortunately, for antiferromagnets, 
i.e. when the magnetic field competes with the spin-spin interaction, this formu- 
lation of the problem leads to a very severe sign problem. In this paper we show 
how this sign problem can be solved completely using a meron-cluster algorithm. 
This method allows us to simulate quantum spin systems in an arbitrary magnetic 
field, in which case the standard loop algorithm fails. Indeed, we show in this paper 
that, at low temperatures or for large values of the magnetic field, the meron-cluster 
algorithm is far more efficient than the loop algorithm. Since the meron-cluster al- 
gorithm has all the advantages of a cluster algorithm, including improved estimators 
for various physical observables, we expect it to be more efficient than the worm 
algorithm. A detailed comparison of the worm algorithm, the stochastic series ex- 
pansion technique and the meron-cluster algorithm would be interesting for future 
studies. 

Antiferromagnetic quantum spin ladders are interesting systems that interpolate 
between 1-dimensional spin chains and 2-dimensional spin systems. The 2-d systems 



may become high-temperature superconductors after doping. As reviewed in flT4fl . 



Haldane's conjecture |15j has been generalized to quantum spin ladders: ladders 
consisting of an odd number of transversely coupled spin 1/2 chains are gapless, while 
ladders consisting of an even number of chains have a gap [16|, [17|, |18|]. Correlations 



in quantum spin ladders were studied numerically in [|T§| , P0| . When the even number 
of coupled chains is increased, the gap decreases exponentially. In the limit of a large 
even number of transversely coupled chains one approaches the continuum limit of 
the 2-d classical 0(3) model |HJ, which can be viewed as a (1 + l)-d relativistic 



quantum field theory. When the quantum spin ladder is placed in an external 
uniform magnetic field B, the corresponding quantum field theory is the (1 + l)-d 
0(3) model with chemical potential \i — B/c (where c is the spin- wave velocity). 
Using Bethe ansatz techniques, this theory has been solved exactly. Translating the 
field theoretic results back into the language of condensed matter physics, we obtain 



2 This was pointed out to one of the authors by H. G. Evertz and M. Troyer some time ago. 
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an analytic expression for the magnetization of a quantum spin ladder consisting 
of a large even number of transversely coupled spin 1/2 chains. The numerical 
results obtained with the meron-cluster algorithm are in agreement with the analytic 
predictions. 

The rest of this paper is organized as follows. Section 2 contains a general dis- 
cussion of the nature of the sign problem and illustrates the basic ideas behind the 
meron-cluster algorithm. In section 3 the path integral representation for quantum 
magnets is derived. Section 4 contains a description of the cluster algorithm for 
quantum spin systems and section 5 introduces the meron concept and discusses 
the solution of the sign problem in detail. In section 6 the physics of antiferromag- 
netic quantum spin ladders is discussed theoretically and the results of numerical 
simulations are presented in section 7. Finally, section 8 contains our conclusions. 



2 Nature of the Sign Problem and the Meron 
Concept 

When antiferromagnetic quantum spin systems in an external magnetic field are 
simulated numerically, a severe sign problem arises. The general nature of this 
problem was discussed in jl]. For the convenience of the reader and to make this 
paper self-contained we reproduce this discussion here. We consider the partition 
function of some quantum system written as a path integral 

Z = ]TSign[ S ]exp(-S[ S ]), (2.1) 

s 

over configurations s with a Boltzmann weight of magnitude exp(— S[s]) and Signfs] = 
±1. Here S[s] is the action of a system with a positive Boltzmann weight and a 
modified partition function 

Z' = £exp(-S[ S ]), (2.2) 

s 

which does not suffer from the sign problem. An observable 0[s] of the original 
system is obtained in a simulation of the modified model as 

(O) = jY,0[s]$w[s)e*p(-S[s]) = ^zv- (2.3) 
l s \sign; 

The average sign in the modified ensemble is given by 

(Sign) = t^E SignN exp(-S[ a ]) = -| = exp(-{3VAf). (2.4) 

The expectation value of the sign is exponentially small in both the volume V 
and the inverse temperature (3 because the difference Af between the free energy 



4 



densities of the original and the modified systems is always positive. Hence the 
original expectation value (O) — although of order one — is calculated as the ratio 
of two exponentially small modified expectation values (OSign) and (Sign) which 
are very small compared to the statistical error and therefore impossible to measure 
in practice. This difficulty is the origin of the sign problem. We can estimate the 
statistical error of the average sign in an ideal simulation of the modified ensemble 
which generates N completely uncorrelated configurations as 

ASign = ^(Sign 2 ) - (Sign) 2 _ exp(pVAf) (o ^ 

(Sign) v^Sign) 

The last equality results from taking the large (3V limit and by using Sign 2 = 1. 
To determine the average sign with sufficient accuracy one needs on the order of 
N = exp(2/?VA/) measurements. For large volumes and small temperatures this is 
infeasible. 

A naive meron-cluster algorithm simply matches any contribution of —1 with 
another contribution of 1 to give 0, such that only a few unmatched contributions 
of 1 remain. Then effectively Sign = 0, 1 and hence Sign 2 = Sign. This step reduces 
the relative statistical error to 

ASign = y/(Sign) - (Sign) 2 _ exy(/3VAf/2) fo ^ 

(Sign) ^(Sign) 

One gains an exponential factor in statistics, but one still needs to generate N' = 
\fN = exp(pVAf) independent configurations in order to accurately determine 
the average sign. Since the probability to measure a contribution of 1 is (Sign), one 
generates exponentially many contributions of before one encounters a contribution 
of 1. Thus this step alone solves only one half of the sign problem and would at 
most allow us to double the space-time volume (3V^\ In a second step involving a 
Metropolis decision, the full meron-cluster algorithm ensures that contributions of 
are suppressed. This step saves another exponential factor in statistics and solves 
the other half of the sign problem. 

A central idea behind our algorithm is to decompose a configuration of quantum 
spins into clusters which can be flipped independently. A cluster flip changes spin 
up to spin down and vice versa. The clusters can be characterized by their effect 
on the sign. Clusters whose flip changes the sign are referred to as merons. Using 
an improved estimator, any configuration that contains merons contributes to the 
average sign, since the flip of a meron-cluster leads to a sign change and hence to an 
explicit cancellation of the contributions of two configurations with the same weight 
but opposite Sign[s] = ±1. Thus only the configurations without merons make 
non-vanishing contributions. If the clusters can be flipped such that one reaches a 

3 The fact that an improved estimator alone cannot solve the sign problem was pointed out to 
one of the authors by H. G. Evertz a long time ago. 
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reference configuration with a positive sign, configurations without merons always 
have Signfs] = 1. This solves one half of the sign problem. 

The other half of the problem is solved using a reweighting method that sup- 
presses multi-meron configurations. In order to perform this essential step, all ob- 
servables must be measured using improved estimators. Fortunately, this is possible 
for a variety of physical quantities. For example, the uniform and staggered mag- 
netization get non-zero contributions only from the zero-meron sector, while the 
associated susceptibilities receive non-zero contributions from the two-meron sector 
as well. Since in this article we will be interested in the uniform magnetization 
only, we can completely eliminate all configurations containing meron-clusters. In 
principle, this restriction can be achieved with a simple reject step. In this way an 
exponentially large portion of the configuration space of the modified model is elim- 
inated. Configurations with merons would contribute in the modified model that 
does not include the sign, but they do not contribute to the observables of the orig- 
inal model with the sign. Eliminating the exponentially large multi-meron sector of 
the configuration space saves an exponential factor in statistics and solves the other 
half of the sign problem. In order to improve the autocorrelations of the algorithm it 
is advantageous not to eliminate multi-meron configurations completely. In practice, 
a Metropolis step is used instead to suppress strongly the incidence of multi-meron 
configurations while still allowing occasional visits of sectors with a small number 
of merons. 



3 Path Integral for Quantum Magnets 

We consider a system of quantum spins 1/2 on a (/-dimensional cubic lattice with 
site label x and with periodic spatial boundary conditions. In particular, we will 
be interested in ladder systems on a 2-d rectangular lattice of size L x V with 
L 3> L'. The spins located at the sites x are described by operators S x with the 
usual commutation relations 

[S l x , — iS xy €ijkS x . (3-1) 

The Hamilton operator 

H = UASlSl + i + SlSlJ + J'S^J - £ B x Sl (3.2) 

x,i x 

couples the spins at the lattice sites x and x + i, where % is a unit-vector in the 
i-direction. The cases J' = J < and J' = J > correspond to the ferro- and 
antiferromagnetic quantum Heisenberg model respectively, while J' = corresponds 
to the quantum XY-model. We consider a magnetic field B x that points in the 1- 
direction and has an arbitrary x-dependence. Later it will be important that the 
magnetic field is perpendicular to the spin quantization axis, which we will choose in 
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the 3-direction. As a specific example we will consider an antiferromagnetic isotropic 
quantum spin ladder (i.e. J' = J > 0) in a uniform magnetic field B x = B. 

To derive a path integral representation of the partition function we decompose 
the Hamilton operator into 2d + 1 terms 

H = H 1 + H 2 + ... + H 2d+1 . (3.3) 

The various terms take the forms 

Hi — 22 h x ^, H i+d = ^ h Xji , H 2 d+i = 51 (3-4) 

x=(xx,x2,...,x d ) x=(xi,X2,...,x d ) x=(xx,x2,...,x d ) 

Xieven xiodd 

with 

h Xii = J(S l X + i + S 2 x S 2 x+l ) + J'S 3 x S 3 x+p h x = -B x S l x . (3.5) 

Note that the individual contributions to a given Hi commute with each other, but 
two different Hi do not commute. Using the Trotter-Suzuki formula we express the 
partition function as 

Z = Tr[exp(-/?iO] = limTr[exp(-/3if 1 ) exp(-/?# 2 )... exp(-^ 2d+1 )] M . (3.6) 

We have introduced (2d + 1)M Euclidean time slices with e = (3/M as the lattice 
spacing in the Euclidean time direction. We insert complete sets of eigenstates | |) 
and | |) with eigenvalues S x = ±1/2 between the factors exp(— (3Hi). The transfer 
matrix is a product of factors 



exp(-eh x . 



( exp(-eJ'/2) \ 

cosh(eJ/2) -sinh(eJ/2) 

-sinh(eJ/2) cosh(eJ/2) 

V exp(-eJ'/2) / 



(3.7) 



(here we have dropped an irrelevant overall prefactor exp(eJ'/4)) and 

, fh] _( cosh(6 J B x /2) sinh(eB x /2) \ 
exp(-eh x ) - ^ sinh{eRx/2) cosHeBx/2) ) ■ I 3 - 8 ) 

The 4x4 matrix of eq. (|3.71) is represented in the basis | ||), | ||), | ||) and | H) 
of two sites x and x + i, and the 2x2 matrix of eq. (p.8|) is in the basis | |) and | |) 
of the site x. 

The partition function is now expressed as a path integral 

Z = ]TSign[ S ]exp(-S[ S ]), (3.9) 
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over configurations of spins s(x,t) =|, | on a (d+ 1 )-dimensional space-time lattice 
of points (x,t). The Boltzmann factor 

exp(-S[s]) = J] ex.p{-S[s(x,t),s(x + i,t),s(x,t + l),s(x + l,t + l)]} 

x = (x 1 ,x 2 ,...,x d ) 

x\even,t=2dp 

x J] exp{-S[s(x,t),s(x + 2,t),s(x,t + l),s(x + 2,t + l)]}... 

x = (x 1 ,x 2 ,...,x d ) 

X2even,t=2dp+l 

x Yi exp{-S[s(x, t),s(x + d, t),s(x, t + l),s(x + d, t + 1)]} 

x = (x 1 ,x 2 ,-..,x d ) 

x d even,t=2dp+d— 1 

x II exp{-S[s(x,t),s(x + i,t),s(x,t + l),s(x + l,t + l)]} 

x = (x 1 ,x 2 ,...,x d ) 

xiodd,t=2dp+d) 

x II exp{-S[s(x,t),s(x + 2,t),s(x,t + i),s(x + 2,t + l)]}... 

x=(x 1 ,x 2 ,...,x d ) 

x 2 odd,t=2dp+d+l 

x ][ exp{—S[s(x,t),s(x + d,t),s(x,t + l),s(x + d,t + l)]} 

x=(x±,x 2 ,...,x d ) 

x d odd,t=2d(p+l)-l 

x II exp{-S[s(x,t),s(x,t + l)]} 

x = (x 1 ,x 2 ,...,x d ) 

t=2d(p+l) 

(3.10) 

(with p = 0,1, M — 1) is a product of space-time plaquette contributions with 

exp(-5[T, T, T, T]) = eM~S[i, i, I, ID = exp(-eJ72), 

exp(-S[T, |, T, I]) = exp(-S[j, t, I, T]) = cosh(eJ/2), 

exp(-5[T, I, I, T]) = exp(-S[j, |, T, I]) = sinh(e| J\/2), (3.11) 

and time-like bond contributions 

exp(-S[T, T]) = exp(-S[j, |]) = cosh(eB x /2), 

exp(-5[T, I]) = exp(-S[l, fl) = sinh(e|^|/2). (3.12) 



The sign of a configuration, Sign[s], also is a product of space-time plaquette 
contributions Sign[s(a;, t), s(x + i,t), s(x, t + 1), s(x + i,t + 1)] with 

Sign[T, T, T, TD = Sign[T, I, T, I]) = Sign[|, |, |, T]) = Sign[|, |, |, |]) = 1 
Sign[T, |, |, T]) = Sign[|, T, T, |]) = -sign(J), (3.13) 

and time-like bond contributions 

Sign[U]) = Sign[UD = 1 

Sign[T, I]) = Sign[|, t]) = sign(^). (3.14) 
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Figure 1 shows two spin configurations in (1 + 1) dimensions. The first configura- 
tion is completely antiferromagnetically ordered and has Sign[s] = I. The second 
configuration contains one interaction plaquette with configuration [J,, |, f, J,] which 
contributes Sign[j, f, f, |] = —1 for J > 0. In addition, there are two time-like in- 
teraction bonds with configurations [j, |] and [T, !]■ When B x > these contribute 
Sign[|, |] = Sign[f , {] = 1, such that the whole configuration has Sign[s] = —1. 

Sign[s] = 1 Sign[s] = — 1 




Figure 1: Two spin configurations in (1 + 1) dimensions. The shaded plaquettes and 
time-like bonds carry the interaction. Filled dots represent spin up and open circles 
represent spin down. For an antiferromagnetic spin chain in a positive magnetic 
field B x > 0, the second configuration has Sign[s] — — 1. 

The central observable of our study is the uniform magnetization 

M* = £5j. (3.15) 

X 

The expectation value of the magnetization (M 1 ) in the direction of the magnetic 
field is non-zero, while (M 2 ) = (M 3 ) = 0. 

It should be noted that the path integral can be simplified for isotropic ferro- 
magnets (J' = J < 0) or antiferromagnets (J' = J > 0) in a uniform magnetic field 
(B x — B). In that case, the magnetic field interaction BJ2 X commutes with the 
rest of the Hamiltonian and can hence be concentrated in a single Euclidean time- 
step. Then one can work with only 2dM+l time-slices instead of (2d+l)M. Instead 
of M magnetic field interaction time-steps with Boltzmann weights cosh(e5/2) and 
sinh(e-B/2) one then has only one such step with Boltzmann weights cosh(/9-B/2) 
and sinh(/?£?/2). In the simulations performed in this paper we use the method with 
(2d+l)M time-slices. 
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4 Cluster Algorithm for the Modified Model with- 
out the Sign Factor 



The meron-cluster algorithm is based on a cluster algorithm for the modified model 
without the sign factor. Quantum spin systems without a sign problem can be 
simulated very efficiently with the loop-cluster algorithm |7], [8], |J. This algorithm 
can be implemented directly in the Euclidean time continuum ||10| , i.e. the Suzuki- 



Trotter discretization is not even necessary. The same is true for the meron-cluster 
algorithm. Here we discuss the algorithm for discrete time. 

The idea behind the algorithm is to decompose a configuration into clusters which 
can be flipped independently. Each lattice site belongs to exactly one cluster. When 
the cluster is flipped, the spins at all the sites on the cluster are changed from up 
to down and vice versa. The decomposition of the lattice into clusters results from 
connecting neighboring sites on each space-time interaction plaquette or time-like 
interaction bond according to probabilistic cluster rules. A set of connected sites 
defines a cluster. In this case the clusters are open or closed strings. The cluster 
rules are constructed so as to obey detailed balance. To show this property we first 
write the space-time plaquette Boltzmann factors as 

exp(— S[s(x, t), s(x + i,t), s(x, t + 1), s(x + i, t + 1)]) = 

C5 s ( x j) t s(x,t+l)S s ( x+ ij),s{x+i,t+l)Ss{x,t)-s{x+i,t) + ^^s(x,t),s(x+i,t+l)"s{x+i,t),s(x,t+l) + 
Ebs(x,t),s{x,t+1)& s^+i^^x+it+l)^ s {x,t),s{x+i,ty (4-1) 

The various 5-functions specify which sites are connected and thus belong to the 
same cluster. It is possible to construct a cluster algorithm satisfying detailed bal- 
ance if the coefficients A,B,...,E are all non- negative. In that case these coefficients 
determine the relative probabilities for different cluster break-ups of an interaction 
plaquette. [] For example, A determines the probability with which sites are con- 
nected with their time-like neighbors, while B and D determine the probabilities 
for connections with space-like and diagonal neighbors. The quantities C and E 
determine the probabilities to put all four spins of a plaquette into the same cluster. 
This is possible for plaquette configurations [T, j, T, j] or [j, f, j, T] with a probability 
proportional to C and for configurations [f, t? 1\ T] or [ii I, j> I] with a probability 
proportional to E. Inserting the expressions from eq. ( |3.11| ) one finds 

expC-StT, T, T, TD = exp(-S[|, I, I, I]) = exp(-eJ'/2) = A + D + E, 
exp(-£[T, 1, T, |]) = exp(-S[|, t, I, TD = cosh(eJ/2) = A + B + C, 
exp(-S[T, 1, 1, T]) = exp(-S[j, T, T, I]) = sinh(e| J|/2) = B + D. (4.2) 

For example, the probability to connect the sites with their time-like neighbors on 
a plaquette with configuration [T, T, T> T] or [ii l> I) I] is A/ (A + D + E), while the 
4 We thank R. Brower for introducing us to the <5-function method. 
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probability to connect them with their diagonal neighbor is D / (A + D + E) . All sites 
on such a plaquette are connected together and are hence put into the same cluster 
with probability E/(A + D + E). Similarly, the probability for connecting space-like 
neighbors on a plaquette with configuration [f, |, |, |] or [J,, f, f, j] is B/(B + £>) 
and the probability for connecting diagonal neighbors is D/(B + D). Similarly, the 
time-like bond Boltzmann factors are expressed as 

exp(-S'[s(x, t), s(x, t + 1)]) = F5 s ( Xj t), s (x,t+i) + G. (4.3) 

The probability to connect spins with their time-like neighbors is F/(F + G). The 
spins remain disconnected with probability G/(F + G). Inserting the expressions 
from eq.( pj,12[ ) one obtains 



exp(-S[T, T]) = exp(-5[|, I]) = coah(eB x /2) = F + G, 

exp(-5[T, I]) = exp(-S[|, |]) = smh(e\B x \/2) = F. (4.4) 

The cluster rules are illustrated in table 1. 

Eqs. Q4.1| , |4.3j ) can be viewed as a representation of the original model as a random 
cluster model. The cluster algorithm operates in two steps. First, a cluster break- 
up is chosen for each space-time interaction plaquette or time-like interaction bond 
according to the above probabilities. This effectively replaces the original Boltzmann 
weight of a configuration with a set of constraints represented by the 5-functions 
associated with the chosen break-up. The constraints imply that the spins in one 
cluster can only be flipped together. Second, every cluster is flipped with probability 
1/2. When a cluster is flipped, the spins on all sites that belong to the cluster 
are flipped from up to down and vice versa. Eqs.( |4.2j , |4.4|) ensure that the cluster 



algorithm obeys detailed balance. To determine A,B,...,E we distinguish three 
cases. For J' > J > we solve eq.(|4.2|) by 



A = exp{-eJ'/2), B = sinh(eJ/2), C = exp(-eJ/2) - exp(-eJ72), D = E = 0. 

(4.5) 

For — J < J' < J we use 

A = i[exp(-eJ72) + exp(-eJ/2)], B = i[exp(eJ/2) - exp(-eJ'/2)], C = 0, 
D = i[exp(-eJ72) - exp(-eJ/2)], E = 0, (4.6) 
and for J' < J < 

A = cosh(eJ/2), B = C = 0, D = sinh(e| J\/2), E = exp(-eJ'/2) - exp(e| J|/2). 

(4.7) 

For example, for the antiferromagnetic quantum Heisenberg model (J' = J > 0) we 
have 

A = exp(-eJ/2), B = sinh(eJ/2), C = D = E = 0. (4.8) 
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weight 


configuration 


break-ups 


exp(-eJ'/2) 


# # 


1 1 X 1 

A D 


□ 

E 


1 


cosh(eJ/2) 








1 


• o 

• o 


11 = 1 

A B 


□ 

C 


sinh(e|J|/2) 


Q • 


B D 


cosh{eB x /2) 


• 




1 

F G 


sinh(e|£y/2) 


i 


■ 

■ 

G 



Table 1: Cluster break-ups of various plaquette and time-like bond configurations 
together with their relative probabilities A, B, ...,G. Filled dots represent spin up, 
open circles represent spin down, and the fat lines are the cluster connections. 



Consequently, on plaquette configurations [t, T, T, T] or [|, |, |, j] one always chooses 
time-like connections between the sites, and for configurations [f , |, |, T] or [l> 1\ t> i] 
one always chooses space-like connections. For configurations [|, I, T> !] or [| 5 T> j> t] 
one chooses time- like connections with probability p = A/(A + B) = 2/[l + exp(eJ)] 
and space-like connections with probability 1— p = B/(A + B). Similarly, eq.( [4.4|) 
yields 

F = exp(-e|£y/2), G = sinh(e| 5,1/2). (4.9) 

The above cluster rules were first used in a simulation of the Heisenberg antifer- 
romagnet [^] in the absence of a magnetic field. In that case there is no sign problem. 
Then the corresponding loop-cluster algorithm is extremely efficient and has almost 
no detectable autocorrelations. When a magnetic field is switched on the situation 
changes. When the magnetic field points in the direction of the spin quantization 
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axis (the 3-direction in our case) there is no sign problem. However, the magnetic 
field then explicitly breaks the Z(2) flip symmetry on which the cluster algorithm is 
based, and the clusters can no longer be flipped with probability 1/2. Instead the 
flip probability is determined by the value of the magnetic field and by the mag- 
netization of the cluster. When the field is strong, flips of magnetized clusters are 
rarely possible and the algorithm becomes very inefficient. To avoid this problem, 
we have chosen the magnetic field to point in the 1-direction, i.e. perpendicular to 
the spin quantization axis. In that case, the cluster flip symmetry is not affected by 
the magnetic field, and the clusters can still be flipped with probability 1/2. How- 
ever, in several cases one now has a sign problem and the cluster algorithm becomes 
extremely inefficient again. Fortunately, using the meron concept the sign problem 
can be eliminated completely and the efficiency of the original cluster algorithm can 
be maintained even in the presence of a magnetic field. 



5 Meron- Clusters and the Sign Problem 

Let us consider the effect of a cluster flip on the sign. As discussed before, the flip of 
a meron-cluster changes Sign[s], while the flip of a non-meron-cluster leaves Sign[s] 
unchanged. An example of a meron-cluster is shown in figure 2 along with the 
same spin configurations as in figure 1. When the meron cluster is flipped the first 
configuration with Sign[s] = 1 turns into the second configuration with Sign[s] = — 1. 
This property of the cluster is independent of the orientation of any other cluster. 
Since flipping all spins leaves Signfs] unchanged, the total number of meron-clusters 
is always even. 

Sign[s] = 1 Signfs] = — 1 




X 

Figure 2: The same spin configurations as in figure 1 together with a meron-cluster 
represented by the fat black line. The other clusters are not shown. Flipping the 
meron-cluster changes one configuration into the other and changes Sign[s]. 
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The meron concept allows us to gain an exponential factor in statistics. Since 
all clusters can be flipped independently with probability 1/2, one can construct an 
improved estimator for (Sign) by averaging analytically over the 2 N ° configurations 
obtained by flipping the Nc clusters in a configuration in all possible ways. For 
configurations that contain merons, the average Signjs] is zero because flipping a 
single meron-cluster leads to a cancellation of contributions ±1. Hence only the 
configurations without merons contribute to (Sign). The probability to have a con- 
figuration without merons is equal to (Sign) and is hence exponentially suppressed 
with the space-time volume. The vast majority of configurations contains merons 
and contributes an exact to (Sign) instead of a statistical average of contributions 
±1. In this way the improved estimator leads to an exponential gain in statistics. 

Now we want to show that the contributions from the zero-meron sector are 
always positive. With no merons in the configuration it is clear that the sign remains 
unchanged under cluster flips, but it is not obvious that one always has Sign[s] = 1. 
To prove this we note that there is no sign problem in the absence of the magnetic 
field. In that case all configurations have Sign[s] = 1 and all clusters are closed 
loops. In the presence of the magnetic field some of the closed loops are cut into 
open string pieces. However, one can always flip the individual pieces such that 
a typical configuration of the model without the magnetic field emerges. These 
reference configurations have Sign[s] = 1. Since any configuration can be turned 
into such a configuration by cluster flips, and since the sign remains unchanged when 
a non- meron-cluster is flipped, all configurations without merons have Signfs] = 1. 
At this point we have solved one half of the sign problem. 

Before we can solve the other half of the problem we must discuss the improved 
estimator for the magnetization. Since the magnetization operator M 1 = J2x &l 
is not diagonal in the basis we have chosen to quantize the spins (which is along 



the 3-direction), it is not obvious how to measure it. However, as discussed in pz 
quantum cluster algorithms naturally allow us to measure non-diagonal operators. 
In the present case, due to time translation invariance, M 1 is equivalent to a sum 
of local spin-flip operators S% inserted at any time slice and appropriately averaged. 
An open string cluster that contains the point (x, t) is naturally divided into two 
parts, one on each side of the point. Flipping each side while keeping the other side 
fixed produces spin configurations that contribute to the spin flip operator Si at 
time t. Closed loop-clusters, on the other hand, do not decompose into two parts 
and thus do not contribute to this operator. If we construct an improved estimator 
for Si inserted at the time t, a configuration containing more than one meron does 
not contribute to (Si). Since the total number of merons is always even, in there 
exists at least one meron that does not contain the point (x, t). Flipping that meron 
changes the sign but not Si and hence results in an exact cancellation of two equal 
and opposite contributions to (Si). Hence we can focus entirely on the zero-meron 
sector. In this sector the point (x, t) can divide an open string cluster into two 
merons or two non-merons. When both parts are non-merons, flipping one part but 
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not the other yields a contribution of 1/2 to (5^). When both parts are merons, 
on the other hand, one gets a contribution of —1/2. Summing over x,t one obtains 
(M 1 ). 

For an antiferromagnet in a uniform magnetic field B x = B, one can find a simple 
expression for the improved estimator for (M 1 ) in terms of a winding number of 
closed loops which result from joining the open string clusters. In that case a non- 
meron cluster has the property that both end points of an open string cluster have 
the same spin type, either up or down. If all the open string clusters are flipped into 
a reference configuration and joined together to form closed loops, one can define 
Wi for each loop to be the temporal winding number of the loop. We assume that 
the cluster is growing in the positive direction on the sites where the loop is broken 
into open strings. Since the loop is always broken up on the same type of spin, the 
above definition of the winding of a given loop is consistent. Further, if a particular 
loop is not composed of open string clusters then W\ = 0. With this definition of 
Wi it is easy to show that 

m = WfioEtM, (5.1) 
2(d7v,o) 

where N is the number of meron-clusters. Hence, as explained above, the mag- 
netization gets non-vanishing contributions only from the zero-meron sector. One 
should note that W\ can be negative. 

Since the magnetization gets non-vanishing contributions only from the zero- 
meron sector, it is unnecessary to generate any configuration that contains meron- 
clusters. This observation is the key to the solution of the other half of the sign 
problem. In fact, one can gain an exponential factor in statistics by restricting the 
simulation to the zero-meron sector, which represents an exponentially small fraction 
of the whole configuration space. This procedure enhances both the numerator and 
the denominator in eq.( |5.1|) by a factor that is exponentially large in the volume, but 
leaves the ratio of the two invariant. In practice, it is advantageous to occasionally 
generate configurations containing merons even though they do not contribute to 
our observable, because this reduces the autocorrelation times. Still, in order to 
solve the sign problem multi-meron configurations must be very much suppressed. 
To suppress such configurations in a controlled way, we introduce an additional 
Boltzmann factor q < 1 for each meron, which determines the relative weights of 
the iV-meron sectors. We visit all plaquette and time-like bond interactions one 
after the other and choose new pair connections between the sites according to the 
above cluster rules. The suppression factor q is used in a Metropolis accept-reject 
step. A newly proposed pair connection that changes the number of merons from 
iV to N' is accepted with probability p = min[l, q N ~ N ]. It should be stressed that 
adding the Boltzmann factor for the reweighting doesn't change the physics, only 
the algorithm. After visiting all plaquette and time-like bond interactions, each 
cluster is flipped with probability 1/2 which completes one update sweep. After 
reweighting, multi-meron configurations are very much suppressed. This completes 
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the second step in our solution of the sign problem. 



6 Antiferromagnetic Spin Ladders in a Uniform 
Magnetic Field 

Antiferromagnetic spin ladders — sets of several transversely coupled quantum spin 
chains — are interesting condensed matter systems which interpolate between single 
1-d spin chains and 2-d quantum antiferromagnets. The ladders are spatially quasi 
1-d systems whose low-energy dynamics are governed by (1 + l)-d quantum field 
theories, which can be solved analytically using Bethe ansatz techniques. This 
allows us to test field theoretical predictions for such models using condensed matter 
experiments as well as numerical simulations. Here we consider spin ladders in an 
external uniform magnetic field B x = B, which corresponds to a chemical potential 
/i = B/c in the corresponding (1 + l)-d quantum field theory. The Bethe ansatz 
solution of this theory yields predictions for the magnetization of the ladder system, 
which we compare with results of numerical simulations. 

Haldane was first to realize the connection between spin chains and (l + l)-d field 



theories ||15|| . He conjectured that a single 1-d antiferromagnetic quantum spin chain 
of length L at inverse temperature (3 is described by a (1 + l)-d 0(3) symmetric 
quantum field theory with an action 

fP r L c 1 i6 

S[e\ = dt dx [—(d x e- d x e+ -^d t e- d t e) + — e ■ {d x e x d t e)]. (6.1) 
jo Jo Ig c z 47r 

Here e(x, t) is a three-component unit- vector field, g is a coupling constant and 
c is the spin-wave velocity. The vacuum angle is given by 9 = 2ttS, where S is 
the value of the spin, hence quantum spin chains with integer spin have 9 = 
and chains with half-integer spin have 9 = 7r. The (1 + l)-d 0(3) model is an 
asymptotically free quantum field theory with a non-perturbatively generated mass- 
gap m oc exp(— 2ir/g 2 ) at 9 = 0. At 9 = it, on the other hand, the mass-gap 
disappears. This naturally explains why chains of half-integer spins are gapless, 
while chains of integer spins have a gap. For large values of the spin S the coupling 
constant is given by 1/g 2 = S/2, so the mass-gap m oc exp(— nS) goes to zero and 
the system approaches a continuum limit. 

Chakravarty, Halperin and Nelson used a (2 + l)-d effective field theory to 
describe the low-energy dynamics of spatially 2-d quantum antiferromagnets |23| . 



Chakravarty has applied this theory to quantum spin ladders with a large even num- 



ber of coupled spin 1/2 chains pl| . We consider spin ladders with the same value 



of the antiferromagnetic coupling along and between the chains. These systems are 
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described by the action 



rP r L r L ' p s 1 
S[e\= dt dx dy -£[d x e- d x e + d y e • d y e + — d t e ■ d t e\, (6.2) 
Jo Jo Jo 2 cr 

where p s is the spin stiffness and c is the spin-wave velocity. The coupled spin chains 
are oriented in the spatial x-direction with a large extent L, while the transverse 
y-direction has a much smaller extent L' <C L. Here we consider spin ladders with 
periodic boundary conditions in the transverse direction, and we limit ourselves to 
an even number of coupled spin 1/2 chains. The effective action for a ladder with 
an odd number of coupled chains would contain an additional topological term. 

When the spin ladder is placed in a uniform external magnetic field B, the field 
couples to a conserved quantity — the uniform magnetization. Hence, on the level 
of the effective theory, the magnetic field plays the role of a chemical potential, i.e. 
it appears as the time-component of an imaginary constant vector potential. As 
a consequence, the ordinary derivative d t e is replaced by the covariant derivative 
d t e + iB x e and the effective action takes the form 

fl3 r L rL' ps 1 ^ 

S[e\ — / dt dx dy — [d x e-d x e+d y e-d v e-\ — - (d t e+iBxe)-(d t e+iBxe)}. (6.3) 
Jo Jo Jo 2 cr 

For a sufficiently large even number of coupled chains [JJ c/p s ) the ladder 
system undergoes dimensional reduction to the (1 + l)-d 0(3) symmetric quantum 
field theory with the action 

S[e\ = dt dx — — \d x e ■ d x e + —(d t e + iB x e) ■ (d t e + iB x e)] 
Jo Jo 2 & 

fPc rL 1 

= / dot dx — [d x e ■ d x e+ (d ct e + ipxe)- (d ct e + ip x e)\. (6.4) 
Jo Jo 2g A 

The effective coupling constant is given by 1/g 2 = p s L'/c and the magnetic field 
appears as a chemical potential of magnitude p = B/c. 

The (1 + l)-d 0(3) model with chemical potential p has been solved exactly 
by Wiegmann |24| using Bethe ansatz techniques. Using this solution, Hasenfratz, 
Maggiore and Niedermayer have derived an exact result for the non-perturbatively 
generated mass-gap m in this model ||25|| . In particular, for p > m and j3 = L = oo 



the Bethe ansatz yields the free energy density 

f(p) - /(O) = ~ f° d9e(9)m cosh9, (6.5) 

Z7T J~6 () 



where the function e(#) satisfies the integral equation 



e(6) - I a d6' (ei _ e ^ +n2 =p-m cosh 6, (6.6) 



70 
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and #o is denned by the boundary condition e(±6 l o) = 0. For // < mno particles can 
be produced and /(//) = /(0). For chemical potentials slightly above this threshold 
the previous equations imply 

M-f(0) = -l — (»-m) 3/2 (6.7) 



and the resulting particle density is given by 

(N) df(n) _ VZm 
L dfi n 



\J \i — m. (6i 



This quantum field theoretic result can be translated back into the language of 
condensed matter physics. It yields an expression for the magnetization density 



(M 1 ) 



L 




(6.9) 



as a function of the magnetic field B. The square root form was derived earlier in 
1 28]. Our result yields an explicit expression for the prefactor as well. 



In the next chapter we will present results of numerical simulations using the 
meron-cluster algorithm. The numerical results are obtained at small but non-zero 
temperatures and for ladders of large but finite lengths, while the above analytic 
expressions are for [3 = L = oo. In particular, the results in the threshold region 
//«m are very sensitive to finite-size and finite-temperature effects. Fortunately, 
these effects can still be understood analytically. In particular, the Bethe ansatz 
solution reveals that for not too large // the system is equivalent to a dilute gas of 
fermions. Its magnetization density is hence given by 



° /L) -^[l + eMM^y + rn'-B))}- 1 . (6.10) 



L L ~y" V L 



Up to this point we have considered the mass-gap m as an unknown constant. 
Based on results from ]33], j2Dfl , Chakravarty |2T| has expressed the mass-gap of 
the ladder system as 

m = ^exp(-27rp s L'/c)[l + + 0{{^-f)\. (6.11) 

ec 4:7rp s L' p s L' 

For the spin 1/2 quantum Heisenberg model with exchange coupling J on a 2-d 
lattice of spacing a, the values of the spin stiffness p s = 0.1800(5) J and the spin-wave 
velocity c = 1.657(2) Ja are known from very precise cluster algorithm simulations 
j27j. As discussed in J27], the above expression for the mass-gap is expected to be 



accurate for a sufficiently large number V / a > 12 of coupled chains. The simulations 
in the present paper are performed at L'/a = 4 and are thus not expected to be 
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described accurately by this formula. Still, for L'/a = 4 the value m = 0.141(2)/a 
is known from the simulations of p0| . 



There is another interesting phenomenon that occurs at rather large values of 
the magnetic field, namely saturation of the magnetization — all spins follow the 
external field and the system becomes completely ferromagnetically ordered. This 
effect cannot be understood in the framework of the above low-energy effective 
theory because it assumes antiferromagnetic order. Still, it is easy to convince 
oneself that there must be a critical magnetic field beyond which saturation occurs. 
This effect is indeed observed in the numerical simulations described below. 



7 Numerical Results 



We have performed numerical simulations with the meron-cluster algorithm for var- 
ious quantum antiferromagnets in a uniform magnetic field B. To demonstrate the 
efficiency of the algorithm, we have compared it with the standard loop-cluster al- 
gorithm. In case of the loop algorithm, the magnetic field points in the direction 
of the spin quantization axis. Hence the Z(2) symmetry that allows the clusters 
to be flipped with probability 1/2 at B — is explicitly broken by the magnetic 
field. In this case some loop-clusters can be flipped only with very small probability 
and the algorithm becomes inefficient for large values of B. In the meron-cluster 
algorithm, on the other hand, the magnetic field is perpendicular to the spin quan- 
tization axis and all clusters can still be flipped with probability 1/2. Now the sign 
problem arises, but it is solved completely by the meron-cluster algorithm. Figure 
3 compares the thermalization behavior of the magnetization of a 2-d Heisenberg 
antiferromagnet on an 8 x 8 lattice at (3 J — 10 with M = 100 for the loop-cluster 
algorithm and the meron-cluster algorithm. At B — 0.75 J, the loop-cluster al- 
gorithm takes about 20000 sweeps until thermal equilibrium is reached, while the 
meron-cluster algorithm equilibrates much faster. At B = J, the loop-cluster algo- 
rithm already needs more than 100000 equilibration sweeps, but the meron-cluster 
algorithm again has no thermalization problem. It should be pointed out that the 
loop-cluster algorithm works reasonably well for small values of B and has been 
used in an interesting numerical study in [^]. However, for larger values of B the 
meron-cluster algorithm is clearly superior. It would be interesting to compare the 



performance of the meron-cluster algorithm with that of the worm algorithm [II, [L2j 
and the stochastic series expansion technique |13| . 



Some of our data for the magnetization density {M x )a/ L of 2-d antiferromagnetic 
spin systems on an L/a x L'/a lattice at inverse temperature (3 are contained in 
table 2. Figure 4 shows the probability to have a certain number of merons in an 
algorithm that samples all meron-sectors without reweighting. This calculation was 
done on an 8 x 8 lattice at (3 J = 1 and M = 100. For small values of B the zero- 
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Figure 3: Thermalization of the magnetization for the loop-cluster algorithm versus 
the meron- cluster algorithm. The first graph is for a magnetic field B = 0.75J and 
the second graph is for B = J. For B = J the loop algorithm takes more than 100000 
sweeps to reach equilibrium while the meron-cluster algorithm has no thermalization 
problem. 

meron sector and hence (Sign) are relatively large, while multi-meron configurations 
are rare. On the other hand, for larger values of the magnetic field the vast majority 
of configurations have a large number of merons and hence (Sign) is exponentially 
small. The meron-cluster algorithm suppresses the multi-meron configurations and 
concentrates on the zero-meron sector which is the only one that contributes to the 
magnetization. 

Figure 5 shows the magnetization density of antiferromagnetic quantum spin 
ladders consisting of four coupled spin 1/2 chains (i.e. L'/a = 4). The cases 
L/a = 20, (3 J = 15, M = 200 and L/a = 40, (3 J = 24, M = 300 have been in- 



20 



L/a 


L'/a 


PJ 


M 


B/J 


(M l )a/L 


8 


8 


10 


100 


0.75 


0.488(6) 


8 


8 


10 
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1.00 


0.690(6) 


20 


4 


15 


200 


0.10 


0.0048(4) 


20 


4 


15 


200 


0.20 


0.0184(4) 


20 


4 


15 


200 


0.30 


0.0452(8) 


20 


4 


15 


200 


0.40 


0.086(2) 


20 
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15 


200 


0.50 


0.120(4) 


20 
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15 


200 


1.00 


0.324(4) 


20 
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200 


2.00 


0.76(2) 
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200 


3.00 


1.280(8) 
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200 


4.00 


1.93(3) 
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4.20 


2.000(8) 


20 


4 


15 


200 


4.40 


2.000(8) 
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24 


300 


0.10 


0.00104(8) 


40 


4 


24 


300 


0.20 


0.0096(6) 


40 


4 


24 


300 


0.30 


0.042(2) 


40 


4 


24 


300 


0.40 


0.085(3) 


40 


4 


24 


300 


0.50 


0.117(7) 


40 


4 


24 


300 


1.00 


0.332(4) 



Table 2: Numerical results for the magnetization density (M 1 )a/L for various spatial 
sizes L, V , inverse temperatures (3, Trotter numbers M and magnetic field values B. 



vestigated. Within our statistical errors these Trotter numbers give results indis- 
tinguishable from the time-continuum limit. In any case, following it would be 
straightforward to implement the meron-cluster algorithm directly in the Euclidean 
time continuum. Using the values m = 0.141(2)/a from |20] and c = 1. 657(2) Ja 
from the finite volume, non-zero temperature expression of eq. QS.10p that was 
derived for small values of B (represented by the dashed curves in figure 5) describes 
the data rather well in that regime. Using the same values for m and c, the infinite 
volume, zero temperature expression for (M 1 )/L that results from solving eqs.( |6.5| ) 
and ( |6.6| ) (represented by the solid curve in figure 5) is in good agreement with the 
numerical data for intermediate values of B. For B/c slightly above the threshold 
m, the infinite volume, zero temperature expression for (M 1 )/L is given by eq.( p.9|) . 
In the threshold region finite size and finite temperature effects are rather large and 
our numerical data do not fall in the applicability range of eq.( |6.9|) . As expected, 
for large values of B the magnetization per spin saturates at 1/2 (represented by 
the dotted curve in figure 5). Note that (M 1 )a/L = 2 for L'/a = 4 in table 2 
corresponds to a value (M 1 )a 2 /(LL') = 1/2 per spin. It should be stressed that 
the comparison of analytic and numerical results does not involve any adjustable 
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Figure 4: The probability of having a certain number of merons for various values 
of the magnetic field B / J = 0.01, 0.05, 0.2 and 0.3. 



parameters. Of course, one should not expect perfect agreement. Given the fact 
that the simulations were performed at L'/a = 4 while the analytic expressions were 
derived in the large L'/a limit, the agreement between theory and numerical data 
is quite remarkable. 



8 Conclusions 



Antiferromagnetic quantum spin systems in a uniform external magnetic field B 
are very difficult to simulate numerically. For example, the standard loop-cluster 
algorithm, which is very efficient for B = 0, suffers from exponential slowing down 
for large values of the field. Here we have traded this slowing down problem for 
a severe sign problem by a change of the spin quantization axis. Meron-cluster 
algorithms provide a general strategy to deal with sign problems and have previously 
led to a complete solution of severe fermion sign problems |], @, |[] . Here we have 
generalized the meron concept to quantum spins in a magnetic field, which again led 
to a complete solution of the sign problem for these bosonic systems. Our numerical 
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Figure 5: Magnetization density {M 1 )/L of quantum spin ladders consisting of 
L'/a = 4 coupled chains as a function of the magnetic field B. The numerical 
data are for two systems: one of size L/a = 20 at inverse temperature (3 J = 15 
(circles) and the other for L/a = 40 at (3J = 24 (diamonds). The solid curve is 
the infinite volume, zero temperature analytic result, while the two dashed curves are 
finite volume, non-zero temperature analytic results for the two simulated systems 
in the small B region. The dotted curve represents saturation of the magnetization 
per spin at 1/2. 

simulations of quantum spin ladders are in agreement with analytic predictions and 
confirm that ladders consisting of an even number of transversely coupled spin 1/2 
chains are described by a (1 + l)-d 0(3) symmetric quantum field theory. 

In conclusion, the meron concept provides us with a powerful algorithmic tool 
- the meron-cluster algorithm — which can lead to a complete solution of severe 
sign problems. In this paper we have demonstrated that this algorithm allows us 
to simulate bosonic quantum spin systems in an arbitrary magnetic field. The 
next challenge is to construct meron-cluster algorithms for systems that show high- 
temperature superconductivity. A first step in this direction was taken in J3|. In that 
paper a Hamiltonian with the symmetries of the Hubbard model but with enhanced 
antiferromagnetic couplings has been investigated. The enhanced antiferromagnetic 
couplings allowed us to solve the fermion sign problem in that model. It turned out 
that after doping the model undergoes phase separation and does not superconduct. 
It remains to be seen if meron-cluster algorithms can also be constructed for models 
that show high-temperature superconductivity. 
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